{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W0D3_LinearAlgebra/student/W0D3_Tutorial3.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Neuromatch Academy: Week 0, Day 3, Bonus Tutorial\n",
    "# Linear Algebra: Discrete Dynamical Neural Circuits\n",
    "\n",
    "__Content creators:__ Name Surname, Name Surname\n",
    "\n",
    "\n",
    "\n",
    "__Content reviewers:__ Name Surname, Name Surname. \n",
    "\n",
    "__Content editors:__ Name Surname, Name Surname.\n",
    "\n",
    "__Production editors:__ Name Surname, Name Surname.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "#Tutorial Objectives\n",
    "\n",
    "In this tutorial, we will start to gain an intuition for how eigenvalues and eigenvectors can be helpful for understanding dynamical systems. We will focus on a discrete dynamical system consisting of two neurons. \n",
    "\n",
    "By the end of the tutorial, you will:\n",
    "\n",
    "* Predict whether the firing rates of interconnected model neurons will explode or decay based on the eigenvalues of the weight matrix.\n",
    "* Apply ideas from previous tutorials (linear combination, basis vectors, etc) to understand a new concept\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:26.431548Z",
     "iopub.status.busy": "2021-06-01T09:18:26.430915Z",
     "iopub.status.idle": "2021-06-01T09:18:27.084160Z",
     "shell.execute_reply": "2021-06-01T09:18:27.083556Z"
    }
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "\n",
    "# Import only the libraries/objects that you use in this tutorial.\n",
    "\n",
    "# If any external library has to be installed, !pip install library --quiet\n",
    "# follow this order: numpy>matplotlib.\n",
    "# import widgets in hidden Figure settings cell\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:27.088585Z",
     "iopub.status.busy": "2021-06-01T09:18:27.088010Z",
     "iopub.status.idle": "2021-06-01T09:18:27.277002Z",
     "shell.execute_reply": "2021-06-01T09:18:27.276445Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Figure settings\n",
    "import ipywidgets as widgets       # interactive display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:27.288973Z",
     "iopub.status.busy": "2021-06-01T09:18:27.287899Z",
     "iopub.status.idle": "2021-06-01T09:18:27.289526Z",
     "shell.execute_reply": "2021-06-01T09:18:27.289911Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Plotting functions\n",
    "\n",
    "def plot_circuit_responses(u, W, eigenstuff = False, xlim='default', ylim='default'):\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10,10))\n",
    "\n",
    "    # Set up axis limits\n",
    "    if xlim =='default':\n",
    "      extreme = np.maximum(np.abs(np.min(u)), np.max(u))\n",
    "      xlim = [- extreme, extreme]\n",
    "    if ylim == 'default':\n",
    "      extreme = np.maximum(np.abs(np.min(u)), np.max(u))\n",
    "      ylim = [- extreme, extreme]\n",
    "\n",
    "    # Set up look\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    cs = plt.rcParams['axes.prop_cycle'].by_key()['color']*10\n",
    "    ax.set_xlim(xlim)\n",
    "    ax.set_ylim(ylim)\n",
    "\n",
    "    # Set up tracking textz\n",
    "    tracker_text = ax.text(.5, .9, \"\", color='w', fontsize=20, verticalalignment='top', horizontalalignment='left', transform=ax.transAxes)\n",
    "\n",
    "    # Plot eigenvectors\n",
    "    if eigenstuff:\n",
    "      eigvals, eigvecs = np.linalg.eig(W)\n",
    "\n",
    "      if np.abs(eigvals[0]) < np.abs(eigvals[1]):\n",
    "        lc1 = 'c'\n",
    "        lc2 = 'g'\n",
    "      else:\n",
    "        lc1 = 'g'\n",
    "        lc2 = 'c'\n",
    "\n",
    "      ax.plot(np.arange(-10000, 10000)*eigvecs[0, 0], np.arange(-10000, 10000)*eigvecs[1, 0],lc1, alpha=.5, label = r'$\\mathbf{v}_1$')\n",
    "      ax.plot(np.arange(-10000, 10000)*eigvecs[0, 1], np.arange(-10000, 10000)*eigvecs[1, 1], lc2, alpha=.5, label = r'$\\mathbf{v}_2$')\n",
    "\n",
    "      ax.legend()\n",
    "\n",
    "    # Set up scatter\n",
    "    cmap = plt.cm.Blues_r\n",
    "    norm = plt.Normalize(vmin=0, vmax=u.shape[1])\n",
    "    scatter = ax.scatter(u[0, :], u[1, :], alpha=1, c = cmap(norm(np.arange(u.shape[1]))))\n",
    "\n",
    "\n",
    "    ax.set(xlabel = 'Neuron 1 Firing Rate', ylabel = 'Neuron 2 Firing Rate', title = 'Neural firing over time')\n",
    "\n",
    "    fig.colorbar(matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap),\n",
    "                ax=ax, label = 'Time step')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:27.299416Z",
     "iopub.status.busy": "2021-06-01T09:18:27.298335Z",
     "iopub.status.idle": "2021-06-01T09:18:27.300035Z",
     "shell.execute_reply": "2021-06-01T09:18:27.300423Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Helper functions\n",
    "\n",
    "def get_eigval_specified_matrix(target_eig):\n",
    "  \"\"\"Generates matrix with specified eigvals\n",
    "\n",
    "    Args:\n",
    "      target_eig (list): list of target eigenvalues, can be real or complex,\n",
    "              should be length 2 unless you desire repeated eigenvalues\n",
    "              with the same eigenvector, in which case length 1\n",
    "\n",
    "    Returns:\n",
    "      ndarray: 2 x 2 matrix with target eigvals\n",
    "\n",
    "  \"\"\"\n",
    "\n",
    "  # Set up two eigenvectors\n",
    "  V = np.array([[1, 1], [-1, 1]]).astype('float')\n",
    "  for i in range(2):\n",
    "    V[:,i] = V[:,i]/np.linalg.norm(V[:,i])\n",
    "\n",
    "  # Get matrix with target eigenvalues\n",
    "  if type(target_eig[0]) == int or type(target_eig[0]) == float:\n",
    "\n",
    "    if len(target_eig) == 2: # distinct eigvecs (not necessarily distinct eigvals)\n",
    "\n",
    "        D = np.diag(target_eig)\n",
    "        A = V @ D @ np.linalg.inv(V)\n",
    "\n",
    "    else: # repeated with same vec\n",
    "      summed = 2*target_eig[0]\n",
    "\n",
    "      a = summed-3\n",
    "      d = 3\n",
    "      bc = target_eig[0]**2 - a*d\n",
    "      factors = [n for n in range(1, bc+ 1) if bc % n == 0]\n",
    "      b = factors[int(np.floor(len(factors)/2))]\n",
    "      c = bc/-b\n",
    "\n",
    "      A = np.array([[a, b], [c, d]])\n",
    "\n",
    "  elif type(target_eig[0]) == complex:\n",
    "\n",
    "      C = [np.real(V[:,0]), np.real(V[:,1])]\n",
    "      B = np.array([[np.real(target_eig[0]), np.imag(target_eig[0])], [-np.imag(target_eig[0]), np.real(target_eig[0])]]).squeeze()\n",
    "      A = C @ B @ np.linalg.inv(C)\n",
    "\n",
    "  return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Section 1: Defining a neural circuit\n",
    "\n",
    "In previous tutorials, we have looked at static models of postsynaptic neurons based on the responses of presynaptic neurons. \n",
    "\n",
    "Let's now introduce the concept of time. We will chop time up into little bins and look at the activity of neurons in each bin. That is, we will work in a **discrete** time framework. For example, if each bin is 1 second long, we will look at the firing rate of each neuron at intervals of a second.\n",
    "\n",
    "\n",
    "Instead of examining pre- and post- synaptic neurons, we will examine at two neurons in one area that are connected. In our model, the activity of neuron 1 at one time bin depends on the activities of both neurons during the previous time bin multiplied by the respective weights from itself and neuron 2. It might seem weird for a neuron to have a weight to itself - this is abstracting away some biological details but basically conveys how much the neural activity depends on its history. (Throughout this course, we'll see lots of neuron models and how some model biological detail more faithfully while others abstract.)\n",
    "\n",
    "We will refer to the activity of neuron i during time bin j as $a_{i, j}$. The weight from neuron x to neuron y will be $w_{y, x}$. With this helpful notation, we can write an equation for the activity of neuron 1 at time bin t:\n",
    "$$a_{1, t} = w_{1, 1}a_{1, t-1} + w_{1, 2}a_{2, t-1} $$\n",
    "\n",
    "And the symmetric model is true of neuron 2:\n",
    "$$a_{2, t} = w_{2, 1}a_{1, t-1} + w_{2, 2}a_{2, t-1} $$\n",
    "\n",
    "This is already a mess of subscript numbers - luckily we can use matrices and vectors once again and our model becomes: \n",
    "\n",
    "$$\\mathbf{a}_{t} = \\mathbf{W}\\mathbf{a}_{t-1} $$\n",
    "where:\n",
    "$$\\mathbf{W} = \\begin{bmatrix} w_{1, 1} & w_{1, 2} \\\\ w_{2, 1} & w_{2, 2} \\end{bmatrix}, \\mathbf{a}_{t} =   \\begin{bmatrix} a_{1, t} \\\\ a_{2, t} \\end{bmatrix}$$\n",
    "\n",
    "It turns out that this is a **discrete dynamical system**. Dynamical systems are concerned with how quantities evolve other time, in this case our neural firing rates. When we model the evolution of quantities over time using a discrete time framework, it is, unsurprisingly, a discrete dynamical system. We will see continuous dynamical systems (where we embrace the full continuity of time) tomorrow and later in the comp neuro course during W2D2: Linear Dynamics.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coding Exercise 1: Implementing the circuit\n",
    "\n",
    "In this exercise, you will implement the function `circuit_implementation`. Given a weight matrix, initial activities at time 0, and a number of time bins to model, this function calculates the neural firing rates at each time bin.\n",
    "\n",
    "We will use initial firing rates of 1 for both neurons:\n",
    "$$\\mathbf{a}_0 = \\begin{bmatrix}\n",
    "1 \\\\\n",
    "1  \\\\\n",
    "\\end{bmatrix}$$\n",
    "and the weight matrix:\n",
    "\n",
    "$$\\mathbf{W} = \\begin{bmatrix} 1 & 0.2 \\\\\n",
    "0.1 & 1 \\\\ \\end{bmatrix}$$\n",
    "\n",
    "We will look at activity over 30 time steps. As before, we will allow our firing rates to be negative, despite this not being possible biologically.\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:27.306288Z",
     "iopub.status.busy": "2021-06-01T09:18:27.305773Z",
     "iopub.status.idle": "2021-06-01T09:18:27.373188Z",
     "shell.execute_reply": "2021-06-01T09:18:27.373684Z"
    }
   },
   "outputs": [],
   "source": [
    "def circuit_implementation(W, u0, T):\n",
    "  \"\"\" Simulate the responses of N neurons over time given their connections\n",
    "\n",
    "  Args:\n",
    "    W (ndarray): weight matrix of synaptic connections, should be N x N\n",
    "    u0 (ndarray): initial condition or input vector, should be N,\n",
    "    T (scalar): number of time steps to run simulation for\n",
    "\n",
    "  Returns:\n",
    "    u (ndarray): the neural responses over time, should be N x T\n",
    "\n",
    "  \"\"\"\n",
    "\n",
    "  # Compute the number of neurons\n",
    "  N = W.shape[0]\n",
    "\n",
    "  # Initialize empty response array and initial condition\n",
    "  u = np.zeros((N, T))\n",
    "  u[:, 0]  = u0\n",
    "\n",
    "  #################################################\n",
    "  ## TODO for students ##\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError(\"Student exercise: Complete circuit_implementation\")\n",
    "  #################################################\n",
    "\n",
    "  # Loop over time steps and compute u(t+1)\n",
    "  for i_t in range(1, T):\n",
    "      u[:, i_t] = ...\n",
    "\n",
    "  return u\n",
    "\n",
    "\n",
    "# Define W, u0, T\n",
    "W = np.array([[1, .2], [.1, 1]])\n",
    "u0 = np.array([1, 1])\n",
    "T = 30\n",
    "\n",
    "# Get neural activities\n",
    "u = circuit_implementation(W, u0, T)\n",
    "\n",
    "# Visualize neural activities\n",
    "plot_circuit_responses(u, W)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:27.379234Z",
     "iopub.status.busy": "2021-06-01T09:18:27.378726Z",
     "iopub.status.idle": "2021-06-01T09:18:27.751489Z",
     "shell.execute_reply": "2021-06-01T09:18:27.751896Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W0D3_LinearAlgebra/solutions/W0D3_Tutorial3_Solution_6fe798cc.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=682 height=702 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W0D3_LinearAlgebra/static/W0D3_Tutorial3_Solution_6fe798cc_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The firing rates of both neurons are exploding to infinity over time. Let's now see what happens with a different weight matrix:\n",
    "\n",
    "\n",
    "$$\\mathbf{W} = \\begin{bmatrix} 0.2 & 0.1 \\\\\n",
    "1 & 0.2 \\\\ \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:27.790561Z",
     "iopub.status.busy": "2021-06-01T09:18:27.775173Z",
     "iopub.status.idle": "2021-06-01T09:18:28.115855Z",
     "shell.execute_reply": "2021-06-01T09:18:28.115432Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to visualize activity over time\n",
    "\n",
    "# Define W, u0, T\n",
    "W = np.array([[.2, .1], [1, .2]])\n",
    "u0 = np.array([1, 1])\n",
    "T = 30\n",
    "\n",
    "# Get neural activities\n",
    "u = circuit_implementation(W, u0, T)\n",
    "\n",
    "# Visualize neural activities\n",
    "with plt.xkcd():\n",
    "  plot_circuit_responses(u, W)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that with this weight matrix, the firing rates are decaying towards zero. It turns out that we could have predicted this by looking at the eigenvalues of the weight matrices, as we'll see in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 2: Understanding dynamics using eigenstuff\n",
    "\n",
    "As we'll see in this section, eigenvectors and eigenvalues are incredibly useful for understanding the evolution of the neural firing rates, and discrete dynamical systems in general.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2.1: Rewriting our circuit equation\n",
    "\n",
    "\n",
    "In our neural circuit, we are modeling the activities at a time step as:\n",
    "$$\\mathbf{a}_{t} = \\mathbf{W}\\mathbf{a}_{t-1} $$\n",
    "\n",
    "Let's start at time  step 1:\n",
    "$$\\mathbf{a}_{1} = \\mathbf{W}\\mathbf{a}_{0} $$\n",
    "\n",
    "And move on to time step 2:\n",
    "$$\\mathbf{a}_{2} = \\mathbf{W}\\mathbf{a}_{1} $$\n",
    "\n",
    "In the above equation, we can subsitute in $\\mathbf{a}_{1} = \\mathbf{W}\\mathbf{a}_{0}$:\n",
    "$$\\mathbf{a}_{2} = \\mathbf{W}\\mathbf{W}\\mathbf{a}_{0} = \\mathbf{W}^2 \\mathbf{a}_{0}$$\n",
    "\n",
    "We can keep doing this with subsequent time steps:\n",
    "$$\\mathbf{a}_{3} = \\mathbf{W}\\mathbf{a}_{2} = \\mathbf{W}\\mathbf{W}^2 \\mathbf{a}_{0} = \\mathbf{W}^3\\mathbf{a}_{0}  $$\n",
    "$$\\mathbf{a}_{4} = \\mathbf{W}\\mathbf{a}_{3} = \\mathbf{W}\\mathbf{W}^3 \\mathbf{a}_{0} = \\mathbf{W}^4\\mathbf{a}_{0}  $$\n",
    "\n",
    "This means that we can write the activity at any point as:\n",
    "$$\\mathbf{a}_{i} = \\mathbf{W}^i\\mathbf{a}_{0}  $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2.2: Initial firing rates along an eigenvector\n",
    "\n",
    "Remember from the last tutorial, that an eigenvector of matrix $\\mathbf{W}$ is a vector that becomes a scalar multiple (eigenvalue) of itself when multiplied by that matrix:\n",
    "\n",
    "$$\\mathbf{W}\\mathbf{v} = \\lambda\\mathbf{v}$$\n",
    "\n",
    "Let's look at what happens if the initial firing rates in our neural circuit lie along that eigenvector, using the same substitution method as in the previous section:\n",
    "$$\\mathbf{a}_{0} = \\mathbf{v} $$\n",
    "$$\\mathbf{a}_{1} = \\mathbf{W}\\mathbf{a}_0 = \\mathbf{W}\\mathbf{v} = \\lambda\\mathbf{v} $$\n",
    "$$\\mathbf{a}_{2} = \\mathbf{W}\\mathbf{a}_1 = \\mathbf{W}\\lambda\\mathbf{v} = \\lambda\\mathbf{W}\\mathbf{v} = \\lambda^2\\mathbf{v}$$\n",
    "$$\\mathbf{a}_{3} = \\mathbf{W}\\mathbf{a}_2 = \\mathbf{W}\\lambda^2\\mathbf{v} = \\lambda^2\\mathbf{W}\\mathbf{v} = \\lambda^3\\mathbf{v}$$\n",
    "$$...$$\n",
    "$$\\mathbf{a}_i = \\lambda^i\\mathbf{v}$$\n",
    "\n",
    "The activities at any time step equal a scalar times the initial activities. In other words, if the initial activities lie along an eigenvector, the activities will only evolve along that eigenvector. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive demo 2.2: Changing the eigenvalue\n",
    "\n",
    "Let's visualize what happens if the initial activities of the neurons lie along an eigenvector and think about how this depends on the eigenvalue.\n",
    "\n",
    "The interactive demo below is the same visualization you saw in Section 1, but now we also plot the eigenvectors $\\mathbf{v}_1$ and $\\mathbf{v}_2$.\n",
    "\n",
    "Questions:\n",
    "1.  What happens if the eigenvalue is large (2)?\n",
    "2.  What happens if you move the eigenvalue from 2 to towards 0? \n",
    "3.  What happens with negative eigenvalues?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:28.136089Z",
     "iopub.status.busy": "2021-06-01T09:18:28.135539Z",
     "iopub.status.idle": "2021-06-01T09:18:28.707412Z",
     "shell.execute_reply": "2021-06-01T09:18:28.706950Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "\n",
    "@widgets.interact(eigenvalue = widgets.FloatSlider(value=0.5, min=-2, max=2, step=0.2))\n",
    "def plot_system(eigenvalue):\n",
    "\n",
    "  # Get weight matrix with specified eigenvalues\n",
    "  W = get_eigval_specified_matrix([eigenvalue, eigenvalue])\n",
    "\n",
    "  # Get initial condition\n",
    "  u0 = np.array([1, 1])\n",
    "\n",
    "  # Get neural activities\n",
    "  u = circuit_implementation(W, u0, 10)\n",
    "\n",
    "  # Visualize neural activities\n",
    "  plot_circuit_responses(u, W, eigenstuff = True, xlim = [-15, 15], ylim = [-15, 15])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:28.713424Z",
     "iopub.status.busy": "2021-06-01T09:18:28.712952Z",
     "iopub.status.idle": "2021-06-01T09:18:28.715149Z",
     "shell.execute_reply": "2021-06-01T09:18:28.714737Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W0D3_LinearAlgebra/solutions/W0D3_Tutorial3_Solution_72e43d55.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2.3: Other initial conditions\n",
    "\n",
    "We now know that if our initial activities (or initial condition) fall on an eigenvector of $\\mathbf{W}$, the activities will evolve along that line, either exploding to infinity if the absolute value of the eigenvalue is above 1 or decaying to the origin it it is below 1. What if our initial condition doesn't fall along the eigenvector though?\n",
    "\n",
    "To understand what will happen, we will use the ideas of basis vectors and linear combinations from Tutorial 1.\n",
    "\n",
    "Let's assume for now that our weight matrix has two distinct eigenvectors ($\\mathbf{v}_1$ and $\\mathbf{v}_2$) with corresponding eigenvalues $\\lambda_1$ and $\\lambda_2$, and that these eigenvectors form a basis for 2D space. That means we can write any vector in 2D space as a linear combination of our eigenvectors, including our initial activity vector:\n",
    "\n",
    "$$\\mathbf{a}_0 = c_1\\mathbf{v}_1 + c_2\\mathbf{v}_2 $$\n",
    "\n",
    "Let's compute the next time step, using our previous strategy of substitution:\n",
    "$$\\begin{align}\n",
    "\\mathbf{a}_1 &= \\mathbf{W}\\mathbf{a}_0\n",
    "\\\\ &= \\mathbf{W}(c_1\\mathbf{v}_1 + c_2\\mathbf{v}_2) \\\\ &= c_1\\mathbf{W}\\mathbf{v}_1 + c_2\\mathbf{W}\\mathbf{v}_2 \\\\ &= c_1\\lambda_1\\mathbf{v}_1 + c_2\\lambda_2\\mathbf{v}_2 \\end{align} $$\n",
    "\n",
    "All activities can be written as:\n",
    "$$\\mathbf{a}_i = c_1\\lambda_1^i\\mathbf{v}_1 + c_2\\lambda_2^i\\mathbf{v}_2 $$\n",
    "\n",
    " We'll see what this means for our system in the next demo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive demo 2.3: Changing both eigenvalues\n",
    "\n",
    "In the demo below, you can now change both eigenvalues and the initial condition (with `a0_1` setting neuron 1 initial activity and `a0_2` setting neuron 2 initial activity). We will only look at positive eigenvalues to keep things a little more simple.\n",
    "\n",
    "Think each of the following questions through based on the equation we just arrived at and then play with the demo to see if you are correct.\n",
    "$$\\mathbf{a}_i = c_1\\lambda_1^i\\mathbf{v}_1 + c_2\\lambda_2^i\\mathbf{v}_2 $$\n",
    "\n",
    "1.  What will happen when both eigenvalues are greater than 1? Does this depend on initial condition?\n",
    "2.  What will happen when both eigenvalues are less than 1?\n",
    "3. Set eigenvalue1 to 2 and eigenvalue2 to 1.2 and try out different initial conditions. What do you see? Why are you seeing this?\n",
    "4. What happens if one eigenvalue is below 1 and the other is above 1?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:28.739120Z",
     "iopub.status.busy": "2021-06-01T09:18:28.738631Z",
     "iopub.status.idle": "2021-06-01T09:18:29.223817Z",
     "shell.execute_reply": "2021-06-01T09:18:29.224196Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to enable the widget\n",
    "\n",
    "@widgets.interact(eigenvalue1 = widgets.FloatSlider(value=0.5, min=0.2, max=2, step=0.2),\n",
    "                  eigenvalue2 = widgets.FloatSlider(value=0.5, min=0.2, max=2, step=0.2),\n",
    "                  a0_1 = widgets.FloatSlider(value=1, min=-5, max=5, step=0.2),\n",
    "                  a0_2 = widgets.FloatSlider(value=2, min=-5, max=5, step=0.2), )\n",
    "def plot_system(eigenvalue1, eigenvalue2, a0_1, a0_2):\n",
    "\n",
    "  # Get initial condition\n",
    "  a0 = np.array([a0_1, a0_2])\n",
    "\n",
    "  # Get weight matrix with specified eigenvalues\n",
    "  W = get_eigval_specified_matrix([eigenvalue1, eigenvalue2])\n",
    "\n",
    "  # Get neural activities\n",
    "  u = circuit_implementation(W, a0, 10)\n",
    "\n",
    "  # Visualize neural activities\n",
    "  plot_circuit_responses(u, W, eigenstuff = True, xlim = [-15, 15], ylim = [-15, 15])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-06-01T09:18:29.229991Z",
     "iopub.status.busy": "2021-06-01T09:18:29.228928Z",
     "iopub.status.idle": "2021-06-01T09:18:29.230542Z",
     "shell.execute_reply": "2021-06-01T09:18:29.230932Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W0D3_LinearAlgebra/solutions/W0D3_Tutorial3_Solution_2993f27c.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2.4: Complex eigenvalues\n",
    "\n",
    "We've been hiding some complexity from you up until now, namely that eigenvalues can be complex. Complex eigenvalues result in a very specific type of dynamics: rotations.\n",
    "\n",
    "We will not delve into the proof or intuition behind this here as you'll encounter complex eigenvalues in dynamical systems in W2D2: Linear Dynamics. \n",
    "\n",
    "Instead, we will simply demonstrate how the nature of the rotations depends on the complex eigenvalues in the animation below. We plot a 3-neuron circuit to better show the rotations. We illustrate each of the following:\n",
    "\n",
    "\n",
    "*   Complex eigenvalues with an absolute value equal to 1 result in a sustained rotation in 3D space.\n",
    "\n",
    "*   Complex eigenvalues with an absolute value below 1 result in a rotation towards the origin.\n",
    "\n",
    "*   Complex eigenvalues with an absolute value above 1 result in a rotation towards the positive/negative infinity.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![complex_eigenvalues.gif]()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "You have seen how we can predict what happens in a discrete dynamical system with an update rule of:\n",
    "$$ \\mathbf{a}_t = \\mathbf{W}\\mathbf{a}_{t-1}$$\n",
    "\n",
    "The most important takeaway is that inspecting eigenvalues and eigenvectors enables you to predict how discrete dybamical systems evolve. Specifically:\n",
    "\n",
    "*  If all eigenvalues are real and have absolute values above 1, the neural activities explode to infinity or negative infinity. \n",
    "\n",
    "*   If all eigenvalues are real and have absolute values above 1, the neural activities decay to 0. \n",
    "\n",
    "*   If all eigenvalues are real and at least one has an absolute value above 1, the neural activities explode to infinity or negative infinity, except for special cases where the initial condition lies along an eigenvector with an eigenvalue whose absolute value is below 1. \n",
    "\n",
    "* If eigenvalues are complex, the neural activities rotate in space and decay or explode depending on the amplitude of the complex eigenvalues.\n",
    "\n",
    "* Even finer details of the trajectories can be predicted by examining the exact relationship of eigenvalues and eigenvectors.\n",
    "\n",
    "Importantly, these ideas extend far beyond our toy neural circuit. Discrete dynamical systems with the same structure of update rule are common. While the exact dependencies on eigenvalues will change, we will see that we can still use eigenvalues/vectors to understand continuous dynamical systems in W2D2: Linear Dynamics. \n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W0D3_Tutorial3",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
